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ABSTRACT 

Massive stars play an important role in many areas of astrophysics, but numerous details 
regarding their formation remain unclear In this paper we present and analyse high resolution 
(R ~ 30 000) near-infrared 2.3 \xm spectra of 20 massive young stellar objects from the RMS 
database, in the largest such study of CO first overtone bandhead emission to date. We fit the 
emission under the assumption it originates from a circumstellar disc in Keplerian rotation. 
We explore three approaches to modelling the physical conditions within the disc - a disc 
heated mainly via irradiation from the central star, a disc heated mainly via viscosity, and a 
disc in which the temperature and density are described analytically. We find that the models 
described by heating mechanisms are inappropriate because they do not provide good fits to 
the CO emission spectra. We therefore restrict our analysis to the analytic model, and obtain 
good fits to all objects that possess sufficiently strong CO emission, suggesting circumstellar 
discs are the source of this emission. On average, the temperature and density structure of the 
discs correspond to geometrically thin discs, spread across a wide range of inclinations. Es- 
sentially all the discs are located within the dust sublimation radius, providing strong evidence 
that the CO emission originates close to the central protostar, on astronomical unit scales. In 
addition, we show that the objects in our sample appear no different to the general population 
of MYSOs in the RMS database, based on their near- and mid-infrared colours. The com- 
bination of observations of a large sample of MYSOs with CO bandhead emission and our 
detailed modelling provide compelling evidence of the presence of small scale gaseous discs 
around such objects, supporting the scenario in which massive stars form via disc accretion. 

Key words: stars; early-type - stars: pre-main sequence - stars: formation - stars: circum- 
stellar matter - accretion, accretion discs 



1 INTRODUCTION 

Massive stars are important from stellar to galactic scales. Their 
high temperature and luminosity results in the injection of large 
amounts of ionizing radiation and kinetic energy into their sur- 
roundings, which shapes the local interstellar medium (ISM) and 
may trigger nearby star formation. They also deposit chemically 
enriched material into the ISM via continuous mass loss during 
their lifetime and in supernova explosions. However, the formation 



* Based on observations made with the ESO Very Large Telescope at the 
Cerro Paranal Observatory under programme ID 079.C-0725 
t E-mail: pyj dilSleeds .ac.uk 



mechanisms of m assive stars (M > SMp) a re poorly understood 
(see the review of IZinnecker & YorkdlioOTh . Given their impor- 
tance in stellar and galactic evolution, it is crucial to understand 
how they form. 

The precursors of massive stars, massive young stellar objects 
(MYSOs ), possess a short Kelv in-Helmholtz contraction timescale 
(lO"*"' yr, iMpttram et al.|| 201 1^ and thus reach the main sequence 
and obtain a high luminosity while still enshrouded in their natal 
cloud material. This high luminosity presents a challenge to the- 
ories of massive star formation, particularly whe n considering a 
sca led up v ersion of low mass star formation (as in lShu et all 19871 . 
see lNorberg & Maeder 2()0 0). There is only a short time for the star 
to accumulate sufficient mass before the protostar ionises the sur- 
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rounding material, which, along with significa nt radiation pressure 
from the protostar, may halt further accretion jLarson & Starrfieldl 
ll97lUKahnlll974IWolfire & Cassinellilll987h . However, recent 3D 
hydrodynamic models indicate that discs circumvent the proposed 
barrier to massive star format ion by facilitating the accretion of 
matter on to the central object dKr umholz eTai]|2009l : lKuiper et alj 
ij). In addition, observations have shown that some ultra- 
compact Hii regi ons have outflows, us ually associated with ongo- 
ing accretion (e.g. lKlaassen et alj201 ih . Thus, continued accretion 
must be possible after the star reaches the main-sequence and be- 
gins ionizing its surroundings. 

Confirming the presence of discs around MYSOs presents a 
considerable observational challenge. Such objects are relatively 
rare, and still embedded in molecular cloud material, making them 
optically invisible. T here have been a handful of detections of discs 
around MYSOs (see IPatel et al.|[20o3; Ijimenez-Serra et al.ll2007l : 



iKraus et alj201(]| : ICarrasco-Gonzalez et alj2012[) but the disc prop 
erties are difficult to determine. Observations at longer wavelengths 
(such as the far infrared and submillimetre) only probe disc proper- 
ties at large distances from the central protostar. Furthermore, very 
few studies can be conducted with sufficient angular resolution to 
probe astronomical unit sized scales, which is necessary to study 
the inner regions of circumstellar discs. The exception is the ob- 
servation of the MYSO G310.0135-l-00.3892 with the Very Large 
Telescope I nterferomete r ( VLT I) and AMBER l lPetrov et al]|2007l) 
reported by iKraus et^ | |20IO|) . which achieved a maximum res- 
olution of approximately 10 au. This provided unique information 
on the geometry of the /f-band continuum emitting material. How- 
ever, this study involved only a single object. Therefore, observa- 
tions of a large sample of MYSOs using an alternative method that 
can probe close to the central protostar are required. 

The inner regions of gaseous discs are difficult to access obser- 
vationally, especially as near-infrared (NIR) interferometric stud- 
ies are l imited to isolated objects which are bright in the near- 
infrared jTatulli et alj|2008l : [Wheelwright et alj|2012bl) . Therefore, 
to study the inner discs of MYSOs, we must employ indirect meth- 
ods. The CO molecule is an ideal tracer of these regions because 
the coupled rotational and vibrational excitation causes a distinc- 
tive emission feature, the CO bandhead, so called because they 
appear in bands in low-resolution spectra. The first overtone v = 
2-0 bandhead emission at 2.3 [im occurs in warm (T = 2500- 
5000 K) and dense (n > lO" cm"^) gas. These are the conditions 
expected in the inner regions of accretion discs. This makes CO 
bandhead emission a valuable tool that allows us to trace these 
regions. In addition, because this feature is the result of transi- 
tions across a range of energy levels (and therefore, temperatures), 
it also allows us to probe the physical properties throughout the 
disc. Previous studies of CO bandhead emission have been suc- 
cessful in fitting spectra of young stars with a range of masses un- 
der t he assumption that the emission originates from a circumstellar 
disc |CarJll989l :'Chandler et al Jll995l:lBik & Thil2004l : lBlum et alj 
OJ; 



l2004lDaviesrt al. 2010: Whe elwright et al.ll20 10l) ."but a study in 



volving a significant number of MYSOs has yet to be performed. 

This has been partly due to a lack of a representative sample 
of MYSOs. Early searches f or MYSOs were conducted using the 
IRAS point source catalogue l lMolinari et al.|[l996l ; ISridharan et alj 
l2002h . This suffered from source confusion due to the large beam 
size (2-5 arcmin at 100 \im) and were biased to isolated objects 
away from the Galactic plane. This issue has been addre ssed by 
the Red MSX Source (RMS) survey jLumsden erai]|2002l) . which 
is an unbiased survey of MYSOs throughout the Galaxy. It is 



drawn from the M SX mid-infrared survey of the Galactic plane 
( lEgan et alj|2003l) . which has a resolution of 18 arcsec, allowing 
detection of sources in previously unresolved regions. An exten- 
sive multi-wavelength campaign has been conducted to identify 
contaminant objects such as ultra-compact Hii regions and plan- 
etary nebula jUrquhart et ai]|2007al . l2009al : iMottram et a"i]|2007l . 
l2010l) . finally yielding approximately 500 candidate MYSOs in 
the database. Kinematic distance estimate s to the MYSOs were 
obtained from molecular line observations jUrquhart et Zll2007bl . 
■2008, 2009b, 20ril. l2012h . while bolometric luminosities have been 
determ ined from fits to the M YSOs' spectral energy distributions 
(SEDs; lMottram et alj201 Ibh . 

In this paper we study a subset of the RMS databasfl We 
utilise our extensi ve low resolution spectroscopic survey of RMS 
sources (see, e.g. IClarke. A. j]|2007l ; ICoope g lin press J) to select 
objects for a high resolution spectroscopic study of CO bandhead 
emission in massive young stellar objects. We detect CO emission 
in 20 MYSOs (and two non-MYSOs), which is compared to kine- 
matic models to assess whether the emission originates in circum- 
stellar discs. Finally, we attempt to determine the properties of the 
CO emitting gas and constrain the accretion rates of these objects. 
Section |2] outlines the observations we have performed while Sec- 
tion [3] details our modelling. Section |4] presents our observations 
and model fits to the spectra, along with an analysis of the best 
fitting parameters. Section [5] discusses our findings and Section |6] 
presents our conclusions from this work. 



2 OBSERVATIONS AND SAMPLE SELECTION 

Table [T| presents the observational parameters of the 20 MYSO, 
and two non-MYSO targets in our study that possessed CO emis- 
sion. The data were taken using the CR IRES near-infrared cryo- 
genic spectrograph dKaeufl etal1l2004) on the Very Large Tele- 
scope (VLT) over three nights in June 2007. A spectral resolution 
of ~ 30 000 was achieved (AA = 0.08 nm at A = 2.3 [im) using a 
slit width of 0.6 arcsec. Standard ABBA nodding along the slit was 
used. The seeing conditions varied from 0.8 to 1.2 arcsec. A sin- 
gle pixel element represents 0.011 nm, while a resolution element 
covers approximately seven pixels. Using a central wavelength of 
/Ic = 2.286 [im, the CO emission spans chips three and four of the 
detector. The CO bandhead peak is located on the third detector 
chip. Telluric standard stars of spectral type A, featureless in the 
wavelength range of interest, were observed between science tar- 
gets using an identical instrumental set-up and at similar airmasses 
to ensure similar sky conditions. 

The data were reduced using the ESO provided CRIRES 
pipeline, via the gasgano data organiser (version 2.2.7). Dark cur- 
rent was removed and detector linearity corrections applied, then 
master flat frames and bad pixel maps were used to correct the 
spectra. The final spectra were extracted using the optimal extrac- 
tion algorithm. Wavelength calibration was p erformed using a cross 
correlation with the raxRAN model catalogue jRothman et alj[l993) 
and OH lines. The standard stars were reduced and extracted in 
the same manner. The final spectra were obtained by division of 
the object spectra by their corresponding standard star to remove 
telluric spectral features. To further remove the effect of bad pix- 
els, the spectra were cleaned using a sigma-clipping algorithm that 
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IRAS 08576-4334 
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Figure 1. Spectrum of the first overtone bandhead emission of IRAS 
08576-4334 (sohd black). The gap in the data is due to the spacing between 
chips tirree (left) and four (right) of the detector. The best fitting model is 
shown above (dashed red) and has been shifted upwards for clarity. The 
theoretical wavelengths of the transitions that make up the bandhead aie 
mai'ked on the abscissa, and have been shifted to account for relative mo- 
tion. The panel shows the line shape of the / = 51-50 transition that was 
adopted for the best fitting model. 



removed any pixels with a value above three times the standard de- 
viation of the pre-bandhead portion of the spectrum. This was de- 
termined to be the maximum amount of cleaning that could be per- 
formed without affecting the appearance of real spectral features. 



2.1 Observational Results 

The spectra, presented in Section |4] exhibit a range of bandhead 
shapes and strengths. Ten objects possess a so-called 'blue shoul- 
der', in which there is prominent emission on the shorter wave- 
length side of the bandhead peak. This is caused by doppler-shifted 
rotational lines, and may be indicative of rotational motion in a 
circumstellar disc, or an outflowing wind ijKraus et al. 20o3). The 
other ten objects show steep rises in the bandhead slope. Several 
objects (for instance, G270.8247-01.il 12) possess a residual tel- 
luric feature in the pre-bandhead portion of the spectrum, but this 
did not affect our analysis. 

The signal to noise ratio of the spectra ranges from approx- 
imately 20-150. The average bolometric luminosity of the oh- 
jects is 5 X 10"* L^, typica l of other objects within the RMS Sur- 
vey jMottram et al.ll201 lah . The objects are generally bright in the 
/f-band, but this is a selection effect based on observational re- 
quirements. 

The objects G332.9457-t-02.3855 and G338.5459+02.1 175 
were originally classified as MYSOs, but since the observation date 
have been found to have too low a luminosity for this to be the case, 
and are likely lower mass young stellar objects. However, since 
they both possess strong CO emission we have included them in 
the final sample and discuss them in Appendix IB] 

FigureEshows the bandhead of IRAS 08576-4334, which is 
a good example of the CO bandhead feature due to the prominent 
emission and the high signal-to-noise ratio. Plotted above the data 
is our final best fitting model (discussed in Section [3j. The inset 
shows the line shape of the J = 51-50 transition that is the result 
of the best fitting model, displaying a double peaked profile due to 
the motion in the disc. 



3 MODELLING THE EMISSION 

To determine whether the CO emission originates from circumstel- 
lar discs and to constrain the physical properties of the emitting 
region, we compare the observations to a model of emission from 
a disc. There are several ways in which the p hysical parameters 
within the disc can be modelled. For example. iBik & Thil ( |2004|) 
assume a simple isothermal disc, ISerthoud et aL 1 20071) assume 
the s urface density and temperature of the di sc decrease as power 
laws. ICard ( Il989h and IChandler et all ( 1 19951) determine the tem- 
perature and density structure by balancing the heating and cool- 
ing mechanisms at work on the surface of a thin alpha disc (see 
IShakura&Sunvaevlll973l) which is isothermal in the vertical di- 
rection. 

To fit the emission using a disc model, we first divide the disc 
into 75 radial rings each with 75 azimuthal cells. Each cell is as- 
signed a temperature and surface density calculated from the rele- 
vant disc description (see Section [3TT] l. We assume Keplerian rota- 
tion for the disc, and the velocity of each cell to the line-of-sight is 
determined assuming the disc is at an inclination i to the vertical. 

Our model of the CO emission is bas ed on [Wheelwright et all 

feOld) , which was in turn based on that of iKraus et alj ilOOd) . and 
briefly described here. The population of the CO rotational levels, 
to a maximum of / = 100, for the v = 2-0 vibrational transition 
is determined in each cell according to the Boltzmann distribution, 
which assumes local thermodynamic equilibrium, and a CO/H2 ra- 
tio of 10"''. The intrinsic line shape of each transition is assumed to 
be Gaussian, with a line width of Av. The intensity of emission is 
then calculated from the equation 



h = By(T)(l-e-'^'). 



(1) 



The emission from each cell is weighted by its solid angle sub- 
tended by the cell on the sky, and wavelength shifted with respect to 
the line-of-sight velocity due to the rotational velocity of the disc. 
The emission from all cells is then summed together, smoothed to 
the instrumental resolution, and shifted in wavelength to account 
for the radial velocity of the object to produce the entire CO band- 
head profile for the disc. 

3.1 Disc models 

We have investigated the use of three different approaches to model 
the disc. Here, we discuss their differences. 

3.1.1 Model A 

Our first model, mo del A, is purely analytic in nature (as in 
iBerthoud et alj[2007h and describes the excitation temperature and 
surface number density as decreasing power laws, 



T(r) = T, 



N(r)-- 



(2) 



(3) 



where Tj and N-, are the excitation temperature and surface density 
at the inner edge of the disc R,, and p and q are the exponents de- 
scribing the temperature and surface density gradient, respectively. 
The optical depth, t, is taken to be the product of the absorption co- 
efficient per CO molecule, and the CO column density. Since we are 
considering a geometrically thin disc, the column density is given 
by the surface number density A'. 
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Table 1. Obsei-ved quantities for each object. The bolometric luminosity LboI. 
database unless otherwise stated. 



2MASS K-hanA magnitude and distance rfy„ are all taken from the RMS 
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3.1.2 Model B 

Our second model, model B, is base d on balanci ng the heating and 
cool ing descriptions of a disc as in ICarJ l ll989l) . later updated by 
[chandler et al. (1995). The disc is assumed to be dominated by 
stellar heating, with a contribution from heating via viscosity, and 
thus the temperature considered is the surface temperature of the 
disc. The disc is assumed to be steady state, and thus has a con- 
stant accretion rate throughout. The disc is under local thermody- 
namic equilibrium. Finally, the disc is isothermal in the vertical 
direction, and is heated via absorption of the stellar radiation field 
and the conversion of gravitational potential energy into thermal 
energy via viscosity. The viscosity within the di sc is described by 
the alpha prescription ISh akura & Sunvaevlll973h . The mass accre- 
tion rate of MYSOs is expected to be high (up to 10"^ Mq yr"' ) and 
the discs are likely turbulent, therefore we adopt a value of o- = 1 .0, 
corresponding to a highly viscous disc environment. 

To obtain the excitation temperature and surface number den- 
sity at a certain radius within the disc, we balance the heating and 
cooling rates at the disc surface. This condition is then iterated 
through with increasing temperatures until fulfilled, and the final 
physical properties are adopted for the particular radius. This is re- 
peated throughout the disc. 

3.1.3 Model C 

Our third model, model C, is again based on balancin g the heating 
and cooling mechanisms of a disc, but as described in lVaidva et alj 



( l2009l) . In this disc, the heating via viscosity is assumed to dom- 
inate the temperature structure, and as such the temperature here 
can be considered the midplane temperature of the disc. The disc is 
steady state and we again adopt a = 1.0. 

We again balance the heating and cooling rates, which pro- 
vides the temperature and number density at a specific radius within 
the disc via iteration. The heating due to irradiation from the star is 
not included in this iteration because we are considering the mid- 
plane temperature of the disc, but is instead added to the temper- 
ature reported from this convergence, yielding the final excitation 
temperature. 

In both model B and C, the mass accretion rate directly sets 
the temperature structure throughout the disc. The inner edge of 
the emission region is set to the radius at which the temperature 
reaches 5000 K (where we assume CO is destroyed), and the outer 
edge of the disc is set to where the temperature drops below 1000 K 
(where we assume CO is no longer ro-vibrationally excited). The 
optical depth is taken to be r = KpH wh ere H is the sca l e heig ht of 
the disc, and k is the opacity taken from [Ferguson et al] ( l2005h . 



3.2 The fitting procedure 

The best fitting model is determined using the downhill simplex 
algorithm, implemented by the amoeba routine of the IDL distribu- 
tion. The input spectra are first continuum subtracted (which is as- 
sumed to be linear given the small range in wavelength), and then 
normalised to the peak of the bandhead emission. Model fits are 



© 2012 RAS, MNRAS QOQ.fTlfn] 



CO bandhead emission ofMYSOs 5 



Table 2. Allowed ranges of parameters for the model fitting procedure, for 
disc models A, B and C. Note that density refers to the surface number 
density. 



Parameter 


Used in 


Range 


Mass accretion rate M 


B, C 


10-'-'' <M < IQ-^'-'Moyr-' 


Inclination i 


A, B, C 


< / < 90 ° 


Intrinsic linewidth Av 


A, B, C 


1 < Av < 30kms"' 


Inner disc radius R, 


A 


1 <Ri < 100 R* 


Inner disc temperature Ti 


A 


1000 < r, < 5000 K 


Inner disc density Ni 


A 


10'2 < M < 10^^ cm-2 


Temperature exponent p 


A 


-5 < p < 


Density exponent q 


A 


-5 < ^ < 



Table 3. Best fitting parameters using disc models A, B and C for the spec- 
trum of G012.9090-00.2607 (W33A). 



Disc Model 


M 

(Moyr-') 


Av 

(kms"') 


i 

n 




A 




21 


37 


1.4 


B 


7.76 X 10-^ 


29 


13 


2.8 


C 


2.14 X 10-5 


29 


27 


2.8 



compared to the data using the reduced chi-squared statistic, x^, 
and the enoi in the data is taken to be the standard deviation of the 
flux in the pre-bandhead portion of the spectra. 

The fitting routine is repeated with six starting positions 
spread across the parameter space to avoid recovering only local 
minima in^J. Table [2] shows the ranges over which we search in 
parameter space. The stellar mass, radius and effective tempera- 
ture are fixed parameters, and are calculated from the bolometric 
luminosity of e ach object using inte rpolation of the main sequence 
relationships in lMartins et al.l l l2005r) , unless otherwise stated in Ta- 
ble[T] The free parameters of the fit are, for model A; the inner disc 
radius, temperature and surface density R,,Ti,Ni, the temperature 
and density exponents p, q, the disc inclination / and the intrinsic 
linewidth Av. For models B and C, the free parameters are: the mass 
accretion rate M, disc inclination ; and intrinsic linewidth Av. 

3.3 A test of the disc models - W33A 

The object G012.9090-00.2607 (hereafter W33A) is a well 
known MYSO tha t has been studie d extensively in recent years 
jde Wit et al] 120101 ; [wheelwright et al. 2012 a), and th i is olfe rs a 
useful check of our models. The work of Ide Wit et alj ilOldi de- 
termined the inclination of the system to be 40° < i < 70° through 
modelling of the SED. We have fitted the spectrum of W33A us- 
ing the three disc models described in Section [3Tl to test which is 
most applicable to the circumstellar environment of MYSOs. Fig- 
ure [2] shows the data and the best fits, and Table [3] shows the best 
fitting parameters for each model. 

Model A clearly gives the best fit to these data, as can be 
seen from the reduced chi-squared statistic. All models reproduce 
the peak of the bandhead and blue side slope well. However, only 
model A accurately reproduces the individual line profiles between 
2.296-2.298 (xm with a sufficiently small intrinsic linewidth. Mod- 
els B and C both have similar chi-squared values, twice that of 
model A. Furthermore, the in clination of model B is too low to 
agree with Ide Wit etakl ilOldi . Model C is more sensitive to high 
accretion rates, which allows higher temperatures in the disc to be 
reached. As can be seen, higher mass accretion rates are reported 
than in model B. The hotter disc means that the CO emission region 
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Figure 2. Best fitting models to G012.9090-00.2607 (W33A) using each 
disc model. Best fitting parameters are shown in Table|3] Model A provides 
the best fit as can be seen from the reduced chi squared statistic. Models B 
and C struggle to reproduce the peak to trough variation seen in the spec- 
trum between 2.296 - 2.298 |xm. All models struggle to fit the initial data 
points on chip four. 

is located further out in the disc, and therefore sufi"ers from less ro- 
tational broadening. This is why model C reports a hi gher inclina- 
tion. H owever, this is still too low to be consistent with lde Wit et alj 
(l2010h . 

Model A provides additional parameters in the temperature 
and density exponents, and their respective values at the inner edge 
of the CO emission region. This allows us to effectively change the 
amount of material within the disc, and provides a better fit to these 
data. While models B & C are both based on physical descriptions 
of accretion discs, they may not contain all of the relevant details. 
For example, these models also assume emission from a flat disc, 
but in reality different disc geometries may need to be considered, 
such as flared discs, or discs with discrete vertical layers. For this 
reason, we chose to adopt disc model A as a basis for the fitting 
routine for the other MYSOs, as it allows freedom to account for 
different emission geometries and is not reliant upon possibly inac- 
curate assumptions regarding the temperature and density structure 
of these discs. 

Adopting disc model A, fits to all objects were obtained as 
described in Section |T2l The errors were calculated by holding all 
but one parameters constant at their best fitting value, and varying 
the selected parameter until the difference in was equal to one. 
For some parameters, the 1-sigma error values were beyond the 
range of allowed parameters for our fitting routine. For all objects, 
we attempted the fitting procedure using data from both chips three 
and four of the detector. If a satisfactory fit was not obtained, we 
limited the wavelength range of fitting to only include data from 
chip three and repeated the fitting procedure. Tests with individual 
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objects showed that extrapolation of a fit using only the wavelength 
range of chip three, on to chip four, produced similar results to a 
fit involving both chips. This can be seen in several objects (for 
instance G296.2654-00.3901 and M8E-IR) where the extrapolated 
fit across chip four is consistent with the location of the rotational 
lines that were not included in the fitting range due to a poor signal- 
to-noise at these wavelengths. 



4 FULL SAMPLE RESULTS 

We are able to fit all objects with spectra that have sufficiently 
strong CO emission, found to be above approximately six times 
the rms noise in the pre-bandhead section of the spectrum. We 
fit 8 MYSOs and the 2 lower mass YSOs across the full chip 
three and four wavelength range, and also use only chip three 
(or a portion of chip three) to fit the remaining 9 MYSOs. 
Thi-ee objects (G293. 8947-00.7825, G332.0939-00.4206 and 
G339.6816-01.2058) were determined to have CO emission that 
was too weak, or a signal-to-noise ratio that was too low, to pro- 
vide a reliable fit, but for completeness we include their spec- 
tra in Appendix iBl The two reclassified low mass YSO objects, 
G332.9457+02.3855 and G338.5459-h02.1175, are excluded from 
our analysis here and discussed in Appendix IB] 

Figure[3]presents the spectra and best fitting models using disc 
model A, to each of the MYSOs with CO emission. Objects where 
only chip three has been used for the fitting routine are indicated 
with a greyed-out chip four region. Table |4] shows the best fitting 
model parameters for each of the objects, with associated 1-sigma 
errors. Errors marked with an asterisk should be considered lower 
limits, as the full 1-sigma error value was beyond our allowed pa- 
rameter range. 

In general, across the 7 objects whose inclinations have been 
constrained before, our best fitting parameters agree with the 
previous results, within error margins. Comments and comparisons 
on an object by object basis are presented in Appendix [A] The 
majority of objects have discs beginning within a few astronomical 
units of the stellar surface, and inner disc temperatures close to 
the dissociation temperature of CO (5000 K). The distribution of 
the inner surface densities of the disc has a geometric average of 
iVi = 5.5 ± 8 X 10-° cm"^. It should be noted that even though CO 
should be di ssociated by stella r U V flux at these smal l dista nces, 
we find as in lSik & Thif ( |2004|) and lWheelwright et al.l feOlOl) that 
the density is sufficient for s elf- shielding to occur (A' > lO'^cm"^, 
Ivan Dishoeck & Blacklll988h . 

Figure|4]shows the distribution of the best fitting inclinations 
and the temperature and density exponents of the MYSOs. While 
our inclinations are consistent with most previously published data 
for objects, IRAS 08576-4334 is not, and we discuss this in de- 
tail in Appendix |A] The distribution of the best fitting inclina- 
tions is essentially consistent with the inclinations being random 
(/ = 60°), as the mean inclination value is ; = 55 ± 25°. The 
temperature gradients are skewed toward higher values, but have 
a mean of p = -0.6 ± 0.5, close to the -0.75 suggested for a 
flat accretion disc. We note that two objects show well constrained 
temperature gra dients of -0.43, which are consistent with flared, 
irradiated discs dChiang & Goldreich|[l997h . The surface density 
gradients are more evenly spread across the parameter space, with 
a mean of ^ = -1.5 ± 1.2. This is consistent with the surface 
density gradient for a flat accretion disc, although it is associated 
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Figure 4. Distributions of inclinations (top), and temperature and surface 
density exponents (p and g, respectively, bottom) determined from the fits 
to the 17 MYSOs with detected CO bandhead emission. 

with a large error. The average intrinsic linewidth of the fits is 
A'y = 20 + 12 km S-' . Apart from two objects (G287.3716-l-00.6444 
and IRAS 17441-2910) the linewidths are approximately ten times 
the thermal linewidths expected for CO between 1000-5000K, sug- 
gesting that the emitting material is dominated by macro-turbulent 
motions or infall. 



5 DISCUSSION 

5.1 Disc sizes - the location of the emission 

We find a range of disc sizes from our results. However, in the ma- 
jority of cases, the inner edge of the CO emission region is within 
a few au. The outer edge of the emission region varies much more, 
giving rise to three very large discs. However, the objects with very 
large discs have temperature exponents that are poorly constrained. 
This allows very shallow gradients (as seen in Figure |4l(, which in 
turn produce large outer disc radii, because this is defined as where 
the excitation temperature drops below 1000 K. Thus, these large 
disc sizes are likely not physical and much smaller discs can be 
produced with an exponent that is still within the uncertainty in the 
best fitting value. 

To compare the location of the CO emission to the circumstel- 
lar disc as a whole, we can calculate the dust sublimation radius, 
Rs, in au for each object: 

^s = 1-1 V&f— ^^^^ — ] ( 'l , (4) 

[IOOOLqI U500K/ 

where 2r is the ratio of absorption efficiencies of the dust, 
and Ts is the temperature at which the dust sublimates 
jMonnier & Millan-GabeJ l2002h . We take Qr = 1 and Ts = 
1500 K. This is then compared to both the inner and outer radii 
of the CO emission disc from our model fits. Figure [5] shows his- 
tograms of the ratios between these quantities. 
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Figure 3. Spectra (black) and model fits (red) to the CO emission of the objects in the sample, using disc model A. Data that are grayed out are not included 
in the fitting procedure because of poor quality, but have been included for completeness. For G035. 1979-00.7427 and G305. 2017+00. 2072, small ranges of 
the spectrum at the edges of the detectors were removed because of excessive noise in the data. 
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Table 4. Quantities used for the model fitting, and results. The stellar mass, radius and effective temperature are derived from bolometric luminosity using 
interpolation from the main sequence relationships described in lMartins et al.l j2005l) , unless otherwise stated. The outer disc radius is defined at the point in 
the disc in which the temperature drops below 1000 K, so no error is presented. If an error value is marked with an asterisk, the change in the value of the 
reduced chi-squared statistic was less than one across the allowed parameter range. 
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Ratio 



Figure 5. Ratio of inner and outer CO disc radii, R[, Ro, compared to the 
dust sublimation radius Rs for each object. The majority of objects have 
inner radius well below the dust sublimation radius, and in some cases the 
entire CO disc is inside the sublimation radius. Most of the discs are within a 
few dust sublimation radii. Note that the final bin in the histogram is uneven 
and extends from 7.5-900. 

As expected, we find that the majority of objects, approxi- 
mately 75 per cent, have CO discs whose inner extent is less than 
the dust sublimation radius. Approximately 30 per cent have CO 
discs with an outer extent below the dust sublimation radius. How- 
ever this percentage is increased if we exclude the artificially large 
discs discussed earlier. Concerning the remaining object, there are 
several cases where the inner and/or outer disc radii are only a 
few times larger than the dust sublimation radius. Our treatment 
of the sublimation is simplistic. Factors such as rapid accretion 
rates, back-warming and non-homogeneous dust grain sizes may 



increase the dust sublimation radius to several times the value we 
predict. Therefore, the best fitting disc properties appear consistent 
with physical expectations. Consequently, we suggest that, in gen- 
eral, the best fitting discs can be associated with gaseous accretion 
discs that are close to the central star. As this is the largest sam- 
ple of MYSOs with CO bandhead emission studied to date, this 
provides strong evidence for the existence of small scale accretion 
discs around these objects. 



5.2 Determining mass accretion rates 

As we discuss previously, the use of a physical disc model that 
directly includes the relevant physics needed to determine a mass 
accretion rate is attractive. However, as these models did not fit 
our data well, we were unable to find accretion rates directly from 
the fits. Therefore, we investigated comparing the radial profiles 
of the temperature and density from our fits to those predicted by 
an accretion disc (model C) for M = 10"3-5_io-7-5 M^yf-i^ an 
effort to determine an estimate for the mass accretion rate of each 
object. We find that, for many of the objects, the mass accretion 
rates obtained from the temperature profile and the density profile 
differed by orders of magnitude. Also, several accretion rates were 
difficult to determine due to very different gradients between our fit 
and the accretion disc model, and the values obtained are therefore 
unreliable. On average, we found that the mass accretion rate from 
the temperature gradients suggested M > lO""* Mq yr"' while the 
surface density distribution suggested M < 10"^ Mq yr"'. 
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Figure 3 - continued Spectra (black) and model fits (red) to the CO emission of the objects in the sample, using disc model A. Data that are grayed out are not 
included in the fitting procedure because of poor quality, but have been included for completeness. The range of fitting for M8E-IR was restricted to obtain a 
good fit, see Appendix |A]for full discussion. For IRAS 17441-2910, the portion of the spectrum beyond 2.306 |xm was not used because an artificial drop in 
flux was caused by the detector 



We conclude that our model fits cannot currently be used to 
determine the mass accretion rates of MYSOs. While this has been 
performed i n previous studies of lower mass young stellar objects 
( jCarJ 19891 : IChandler et al.lll995l) . the higher resolution of our data 
demonstrates there are many rotational lines that cannot be fitted 
using these simple disc models. These models may be inadequate 
in accurately representing the physical situation that gives rise to 
CO bandhead emission from discs with high accretion rates. Also, 
the study of only a single bandhead may not offer sufficient infor- 
mation to reliably determine the physical properties of discs and 
thus cannot accurately constrain mass accretion rates in this way. 
Observations covering a larger range in wavelength, for instance 



with VLT/XSHOOTER, would allow fitting of additional band- 
heads that probe different temperature regimes, at the expense of 
the high spectral resolution employed in this paper. 



5.3 Are these typical MYSOs? 

In this section we consider whether the MYSOs with CO emission 
differ from those without. An important question since the detec- 
tion rate of CO emission in the spectra ofMYSOs is approximately 
25 per cent. 

MYSOs are typically red objects with a featureless continuum 
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in the optical and NIR, so we cannot accurately constrain their stel- 
lar properties such as effective temperature and radius. Therefore, 
we investigate whether the parameters of the best fitting models dis- 
play a pattern which may explain why only some MYSOs exhibit 
this emission. In particular, it is conceivable that for CO emission 
to be observed, the inclination of the disc is required to be close 
to face on. However, the inclinations of the objects are spread rela- 
tively evenly between and 90°, with a slight preference for higher 
inclinations, suggesting this is not the case. In addition, we note that 
the average bolometric luminosity of our sample of objects with CO 
emission is 5 x 10"* Lq, which is typical of objects within the RMS 
survey. Therefore, we conclude that the properties of the MYSOs 
that exhibit CO bandhead emission do not indicate what specific 
geometry and/or conditions result are required for the presence of 
this emission. 

Consequently, it is not certain why the presence of CO band- 
head emission is not ubiquitous in the spectra of MYSOs. It is pos- 
sible that the objects with CO emission represent a different evolu- 
tionary stage of MYSOs than those without. It is difficult to test this 
hypothesis. As an initial test, we investigate the infrared colours of 
the observed MYSOs, which are likely affected by key factors such 
as circumstellar geometries on astronomical unit scales, inclination 
and envelope mass/infall rate. To determine whether the MYSOs 
with CO bandhead emission appear representative of MYSOs in 
general or are a specific subset of MYSOs, we compared their 
NIR colours with those of approximately 70 objects from the RMS 
database, which we show in the upper panel of Figure [6] To en- 
sure a valid comparison with our objects, the control sample was 
selected to have a high luminosity (L > IO'^Lq) and be bright in 
the A'-band {K < 10 magnitudes). We compared the J - K and 
H - K colours of the two samples with the Kolmogorov-Smimov 
(KS) test and found that the hypothesis that the NIR colours of the 
MYSOs with CO bandhead emission are drawn from the distribu- 
tion of NIR colours exhibited by the RMS population cannot be 
discounted with any significance. The objects observed to have CO 
bandhead emission are generally bright in the A^-band. However, 
this is a selection effect and does not imply that only objects with 
such emission are intrinsically brighter in the A'-band. 

We also examined the mid-infrared (MIR) colours of the sam- 
ple. Specifically, we compared the MSX {Fji^im - Fiiim), (Piipm ~ 
Fnijm} and (F2i^m - ^i4;jm) colours of the objects with CO bandhead 
emission to the RMS population mentioned earlier. The lower panel 
of Figure|6]shows the cumulative distribution of the {Fii^m - Fsfim) 
colours. Using the KS test, we found that the hypothesis that the 
MIR colours of the objects with CO emission are drawn from 
the total RMS distribution of MIR colours cannot be discounted 
with any significance. Therefore, the objects with and without CO 
emission appear no different in terms of both their NIR and MIR 
colours, suggesting these objects are representative of the RMS 
population as a whole. 

It is therefore unclear from this analysis why only some 
MYSOs possess CO emission. It has been predi cted that models 
of cir cumstellar discs of MYSOs can be unstable jRrumholz et alj 
and that the accretion rate in these discs is not constant 
jKuiper et al. I I20T ij). Unstable discs may disrupt the circumstellar 
material and lead to physical conditions in which CO ro- vibrational 
emission no longer occurs. If the disc accretion rate is high, this 
will increase the mid-plane temperature of the disc and move the 
CO emission region further away from the central protostar, possi- 
bly into the surrounding envelope, which may cause the emission to 




Figure 6. Cumulative frequency of the NIR J - K and normalised MIR 
(Fiifjm ~ f^8(/m) colours of the MYSOs with CO emission and the RMS 
population. The dashed vertical lines indicate the largest deviation between 
the distributions. Similar distributions were found for the H - K, (F2ipm ~ 
Fn^m) and (F^i^m - F'wpm) colours. 

cease. However, without a method to determine the accretion rate 
of these objects, this hypothesis is difficult to test. 

5.4 Discs and massive star formation 

The physical processes occurring during massive star formation 
are still not well understood. It is thought accretion proceeds 
through circumstellar discs at high ra tes, which are believed to af- 
fect the state of the ce ntral protostar ( iHosokawa & Omukaill2009l : 
iHosokawa et alj I2OIOI) . Furthermore, recent simulations indicate 
that high accretion rates result in massive discs and that gravita- 
tional torques in such self-gravitating di scs provide a mech anism to 
transport angular momentum outwards jKuiper et alj201 ll) . There- 
fore, it is becoming apparent that accretion discs play a central role 
in the formation of massive stars. However, it is difficult to confirm 
this observationally. In most cases, the large distances to MYSOs, 
and their embedded nature, prevent detections of astronomical unit 
scale circumstella r discs (the exception being IRAS 13481, see 
iKraus et al]|20IOl) . In some cases, high infall rates have been de- 
tected towards massive star forming regions (see e.g. Beuther et alj 
2012; Wvrow.ski et al. 2012; Oiu et al. 2012; Herpin et al. 201^ 
However, these observations probe scales of approximately 1000 au 
and larger. Therefore, it is difficult to establish that this material 
will be accreted by a single object. Consequently, our observational 
overview of massive star formation is still incomplete. 

With this in mind, we note that our sample constitutes the 
largest sample of MYSO s with CO bandhead emis sion studied to 
date (several times that o f lWheelwright et alj|2010l) . Therefore, by 
confronting the observations with a kinematic model, we can con- 
duct an extensive investigation into the circumstellar environment 
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of these objects. We demonstrate that all the observed bandhead 
profiles can be successfully fitted with a model of a circumstellar 
disc in Keplerian rotation. In addition, we have shown that essen- 
tially all these models can be associated with gaseous discs interior 
to the dust sublimation radii of these objects. Finally, we demon- 
strate that the objects in question appear no different to the ensem- 
ble of objects in the RMS catalogue, the largest most complete cat- 
alogue ofMYSOs to date. Therefore, the fact that the CO bandhead 
emission of all the objects observed can be fitted with a disc model 
supports the scenario in which all MYSOs are surrounded by accre- 
tion discs. With our limited wavelength coverage, we were not able 
to establish precise constraints on the accretion rates. Nonetheless, 
we find that the disc properties are consistent with high accretion 
rates. Therefore, these observations are entirely consistent with the 
hypothesis that MYSOs are massive stars in the process of forming 
that accrete matter through circumstellar discs. 



6 CONCLUSIONS 

In this paper we have presented the near-infrared spectra of 20 mas- 
sive young stellar objects that possess CO first overtone bandhead 
emission, the largest sample studied in this way to date. We have 
fit the spectra of these objects with a model of emission originat- 
ing in a circumstellar disc in Keplerian rotation. We tested three 
approaches to describing the properties of such circumstellar discs 
and found that the spectra were best fit by using an analytic ap- 
proach to describe the temperature and density within the disc. Our 
main findings are: 

(i) All spectra are well fit by a model of a Keplerian rotation 
disc. 

(ii) The best fitting disc parameters are consistent with previ- 
ously published information. The inclinations are spread across a 
wide rage of angles. The best fitting temperature and density expo- 
nents are, on average, consistent with flat circumstellar discs (sub- 
ject to some scatter), a handful of objects have exponents consistent 
with flared, irradiated discs. 

(iii) Essentially all the best fitting discs are located close to the 
dust sublimation radius, which is consistent with the existence of 
small scale gaseous accretion discs around these objects. 

We found that the mass accretion rates of the objects are not 
easily determined from examination of the disc structure as traced 
by a single CO bandhead. The analysis of further bandheads may 
allow these accretion rates to be determined. Our high spectral res- 
olution observations of the CO emission of a substantial number of 
MYSOs and our modelling of this emission are entirely consistent 
with the scenario in which MYSOs are surrounded by small scale 
(< 100 au) accretion discs. The objects observed constitute a large 
sample of circumstellar discs around MYSOs and provide promis- 
ing targets for future observations and inspiration for detailed mod- 
elling of the accretion environment of such objects. 
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Al G033.3891+00.1989 

The work of [Wheelwright et al.l 1 I2OIOI) fitted the CO bandhead 
emission of this object using fixed power law relations for the disc. 
They found a disc from 0.24-2.0 au, at an inclination of 18°, with 
an intrinsic linewidth of 20 km s"' (which was not well constrained) 
and a CO number density of 9 x 10^' cm"-. The inclination, number 
density and linewidth agree with our best fitting model within er- 
ror ranges, but the size of our disc is much larger, due to a shallow 
temperature gradient. 

A2 G287.3716-00.6444 

IWheelwright et all ( I2OI0I) do not find a satisfactory fit to the CO 
bandhead of this object assuming a circumstellar disc (with fixed 
temperature and density profiles), nor an isothermal non-rotating 
body of CO. They discuss other possible origins for the emission, 
including a disc with an outer bulge, a dense neutral wind, a shock, 
or a disc in which the receding side is much brighter than the ap- 
proaching side. We note that our temperature exponent is close 
to -0.43, which would be consistent with a flared, irradiated disc 
jChiang & Goldreicd [19971) which may act in the same way as a 
disc with an outer bulge. 

A3 G308.9176-h00.1231 (AFGL 4176) 

[Wheelwright et al.l 1 I2OI0I) find a disc from 1-8 au at an inclination 
of 30°. The size of the disc agrees well with our results, but we find 
a higher inclination of 67°, with 30° at the lower limit of our er- 
ror range. Our linewidth of 1 2.6 km s"' agrees with their value of 
14km s~', however our inner density is one order of magnitude 
lower. ISolev et al.l l l2012l) find their observations described well by 
a large circumstellar disc at an inclination of 60°, agreeing with our 
best fitting model, and consistent with the prominent blue shoulder 
in our data. 



A4 G310.0135+00.3892 

G3 10.0135+00.3892 (IR AS 13481-6124) ha s previously been ob- 
served in the A^-band bv lKraus et all ( 1201 Ol) who report an elon- 
gated structure that is consistent with a disc viewed at a mo derate 
inclination of approximately 45°. [Wheelwright et al.l l l2012ah fitted 
the SED using a model with an inclination of 32°. Our relatively 
high inclination of 67° is not well constrained, due to the poor qual- 
ity of the data, but agrees with these values within the error range. 

iKraus et al] 1 I2OIOI) find a temperature gradient of p ~ -0.4, 
which they sug gest is consistent with a flar ed, irradiated disc based 
on the work of Chiang & GoldreichI ( Il997l) . We find a temperature 
gradient of p = -0.43, which is consistent with this hypothesis. 
Our inner disc temperature of 3800 K is warmer than the v alue of 
approximately 1500-2000 K assumed in lKrausetal.l ( l2010h . but is 
consistent as we are concerned with a gaseous disc as opposed to a 
dust disc. We find a smaller inner radius (2.8 au) for our disc than 
their study (9.5 au). 



APPENDIX A: NOTES ON INDIVIDUAL OBJECTS & 
COMPARISON OF RESULTS WITH PREVIOUS STUDIES 

Here we discuss and compare, on an object by object basis, our 
findings with that of previous studies on a selection of our sample 
where data was available. 



A5 G332.9868-00.4871 

[wheelwright et all ( l2012al) determine an inclination of 15° to 
G332.9868-00.4871, which is far from our reported value of 78° 
even considering the large error. However, we note that our data 
have a relatively low signal-to-noise ratio, and no rotational lines 
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in the fourth chip can be observed, thus our fit is likely not the best 
fit for the object. 

A6 G347.0775-00.3927 

IWheelwright et aljfeoioh find a similar sized disc to our best fitting 
model, from 0.5^ au but with a lower inclination of 30°. However, 
our inclination of 84° is not well defined, and this does lie within 
our lower error bound. It should also be noted that the bolometric 
luminosity from the RMS database used in their study has since 
been revised to a lower value, which we use here. 

A7 IRAS 08576-4334 

IRAS 08576-43 34 has been extensively studied in recent years. 
iBik & Thil Jaooj) model the CO emission using an isothermal disc 
from 0.2-3.6 au, with an inclination angle of 27°, a CO number 
density of 3.9 x 10^' cm~^ and an excitation temperature of 1600 K. 
IWheelwright et alj feoioh used a similar method, but utilise a disc 
model with fixed power laws, and show the data to be well fit with 
a 0.09-0.78 au disc, at an inclination of 18°, with a CO number 
density of 7.9 x 10^' cm"^ and an intrinsic linewidth of 20 km s ' . 

lEllerbroek et al. I( |20T ij) find double peaked He i and Fe i emis- 
sion lines, with a separation of 60-100 km s"', which they con- 
clude must originate from a circumstellar disc. They suggest IRAS 
08576-4334 is likely an intermediate mass YSO with a mass ac- 
cretion rate of 10~*-I0~' Mg yr~', obtained from the determination 
of the outflow mass loss rate. 

In contrast to lBik & ThU ( l2004h and lWheelwright et alj ilOldi . 
we find our observations are best fitted with larger disc from 
0.6-6.5 au, at an inclination of 65°. We note that data of 
[Wheelwright et alj ilOld) does not have sufficient wavelength 
range to observe double-peaked emis sion lines beyond 2.297 yun, 
and that the lower resolution of the iBik & Thil ( l2004b data may 
mask the presence of these features, especially as they are not well 
defined and asymmetric in our data. Our best fitting model pos- 
sesses a double peak width of approximat ely 50kms~', similar t o 
that of the He i and metal emission lines in lEllerbroek et"ZI ( l20Ill) . 

We note that the rotationally broadened lines from 2.306- 
2.309 vim show asymmetry, with a depletion on the blue side. These 
transitions likely correspond to the cooler material further out in the 
disc, and as such may be evidence for asymmetry in the disc, which 
cannot be fitted with our axisymmetric disc model. 

A8 IRAS 16164-5046 

iBik & Thil t004) model the CO emission of IRAS 16164-5046 
using an isothermal disc from 3. 1-3.2 au (with large errors) at an 
inclination of 30° , with an excitation temperature of 4480 K and a 
number density of 4x lO-^cm"'. The extent of the disc is consistent 
with our results within errors, however we find a best fit that is 
closer to the central star, which would account for our larger inner 
number density. Our disc inclination of 53° is higher. 

iBik et al.] ( |2006|) find CO emission and Pfund series emission 
with a width comparable to that of the CO emission, suggesting 
a common kinematic origin. However, due to the different con- 
ditions required for both emission, they suggest the CO emission 
comes from the midplane of a disc, while the Pf-lines comes from 
the ionised upper layers. They also note CO absorption from 2.33- 
2.35 \im, indicative of cold, foreground molecular gas. 



A9 M8E-IR 

IWheelwright et al.l ilOvh were able to fit the emission of this ob- 
ject with a disc from 0.3-2.6 au at an inclination of 16°, with an 
inner number den sity of I x 10^' cm"^ and a linewidth of 7 km s"' . 
iLinz et al. I ( l2009h find an inclination of 18.5° and a density expo- 
nent of ^ = -2.05 using a = 0.013. 

Our fits to M8E-IR did not satisfactorily converge using many 
starting positions across the initial parameter space. The spectrum 
was highly reddened, and the fit suffered from several minima with 
similar values. The only change applied between these fits was 
the level of continuum subtraction applied. There were two issues 
with the data. Firstly, the level of the chip 4 flux seemed too high 
to be reproduced by the model, meaning that even solutions that 
reproduced the rotational line structure were assigned poor^^ val- 
ues. Secondly, the model was unable to reproduce the relatively 
narrow structures on the red side of the bandhead edge, which may 
be noise. Because these features were across a larger range of wave- 
length than the three broader lines at the edge of chip 3, the fitting 
routine assigned them a higher weight which resulted in poor fits to 
the bandhead slope and the broader lines. 

To address this, we restricted the range of fitting to exclude 
this region (as can be seen in Figure |4ll, and increased the allowed 
upper linewidth to 60kms"', which produced a better fit to the 
d ata, and produced similar best fitting parameters similar to those 
in lLinz et al] ( l2009h . 



APPENDIX B: NON-MYSO SOURCES & OBJECTS THAT 
WERE NOT FITTED 

The objects G332.9457+02.3855 and G338.5459+02.1 175 were 
originally thought to be MYSOs at the time of observing, but subse- 
quent determination of their bolometric luminosity has shown the 
objects too faint for this to be the case, and they are likely lower 
mass young stellar objects. We could not base our estimation of 
their stellar parameters on the main sequence relationships, so we 
have estimated their stellar parameters at the values shown in Table 
151 which lie in the range of typical T Tauri star values. Changes to 
these parameters within these ranges had little effect on the final 
fits. The fits are presented in Figure IBTI 

Magnetically channelled accretion funnels have been sug- 
geste d as a possibl e source for CO bandhead emission in T Tauri 
stars ( lMartin|[l997l) . We obtain good fits to the data of two young 
stellar objects using a simple disc model. The sizes of the discs are 
beyond the typical co-rotation radii for these objects. In addition, 
the intrinsic linewidths of both objects are similar to those of the 
massive YSOs. This suggests that the emission originates from cir- 
cumstellar discs regardless of the mass of the central YSO, and not 
from accretion funnels. 

Finally, we note that three objects in the sample had CO emis- 
sion that was too weak, or a signal-to-noise ratio that was too low 
for an accurate fit to be obtained. For completeness, we include 
their spectra in Figure |B2| 

This paper has been typeset from a TgX/ KTgX file prepared by the 
author. 
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Figure Bl. Objects that were later confirmed to be non-MYSO sources, but 
possess strong CO emission and were fitted with our disc model. 
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Figure B2. Objects that possessed CO emission that was too weak for an 
accurate fit to be obtained. 
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